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We numerically study the vortex-vortex interaction in multi-component homogeneous Bose- 
Einstein condensates within the realm of the Gross-Pitaevskii theory. We provide strong evidences 
that pairwise vortex interaction captures the underlying mechanisms which determine the geometric 
configuration of the vortices, such as different lattices in many-vortex states, as well as the bound 
vortex states with two (dimer) or three (trimer) vortices. Specifically, we discuss and apply our the¬ 
oretical approach to investigate intra- and inter-component vortex-vortex interactions in two- and 
three-component Bose-Einstein condensates, thereby shedding light on the formation of the exotic 
vortex configurations. These results correlate with current experimental efforts in multi-component 
Bose-Einstein condensates, and the understanding of the role of vortex interactions in multiband 
superconductors. 
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I. INTRODUCTION 

The realization of Bose-Einstein condensates (BECs) 
has brought about a suitable research environment for 
investigating general properties of superfluidity and su¬ 
perconductivity with a high degree of control and ver¬ 
satility m- This has led vortex states and their dy¬ 
namics, key concepts for both superfluidity and super¬ 
conductivity, to the rank of some of the most investi¬ 
gated topics in low temperature physics [IHB]. Remark¬ 
ably, since the hrst observation of vortices [Hi and, 
subsequently, the formation of highly-ordered vortex lat¬ 
tices in BECs [H], much theoretical and experimental ef¬ 
fort has been made to push further the understanding 
of these systems. Specially, recent advances in experi¬ 
mental techniques have led to the realization of conden¬ 
sates with several different types of particle-particle in¬ 
teraction, making this research field even broader. In 
the spotlight, theoretical results on dipolar mm and 
spin-orbit-coupled [T^] BECs exhibit rich ground-state 
phases, with bubbles, density stripes and various vor¬ 
tex lattice geometries. Besides the already mentioned 
possibilities, further remarkable physical phenomena to 
be addressed with ultracold atomic systems include, for 
example, quantum-fluid turbulence in BECs m, mul- 
ticharged vortices [14] and vortex-antivortex lattices in 
superfluid fermionic gases [15]. 

With the recent prospects of producing condensates 
with a large number of components and different types 
of interaction [Min], vortex-states in multi-component 
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condensates become even more important. Indeed, recent 
theoretical results on vortex lattice conformations have 
shown that multi-component Bose-Einstein condensates 
are far from being just a trivial extension of the single 
component case. In fact, adding a component brings a 
diversity of possible configurations never found in a one- 
component system, such as amorphous conformations, 
square lattices and bound states, including overlapped 
vortices, vortex dimers and molecules [THH13- Some of 
these conformations suggest that the interaction between 
vortices can be nonmonotonic with respect to the inter¬ 
vortex distance. Unfortunately, to date, the interaction 
between vortices in the simplest multi-component BEC 
is known only in specihc limits of scale, by either assum¬ 
ing inter-vortex distances much greater than the heal¬ 
ing length |23j or considering the interaction energy near 
the vortex peak in the Thomas-Fermi regime [23]. As it 
turns out, the asymptotic behavior does not account for 
all conformations found [I9l - l22| . Furthermore, the gener¬ 
alization of these analytical approaches to more complex 
cases with more components or even with a different kind 
of inter-component particle-particle interaction seems to 
be highly non-trivial. 

Within this context, we investigate in the present pa¬ 
per the origin of these unusual quantized vortex states 
by focusing on the pairwise vortex-vortex interaction 
[2H1 I2H] . We consider homogeneous BECs, which pos¬ 
sess translation invariance. On the one hand, this is im¬ 
portant in itself, since the key properties of the exper¬ 
imentally more relevant harmonically trapped BECs in 
the large particle number regime bear close resemblance 
to those of their homogeneous counterparts m- On the 
other hand, the recent achievement of a condensate in 
a uniform potential [28] enables the experimental veri¬ 
fication of our predictions. The experimental progress 
has already led to the detection of extended phase coher¬ 
ence in a uniform quasi-two-dimensional Bose gas [29] . 
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In our approach, constraints are imposed on the Gross- 
Pitaevskii (GP) formalism, allowing for fixing the vor¬ 
tices in desired positions for further analysis. In other 
words, vortices no longer arise naturally from the GP 
equation itself, but are instead placed manually in the 
position of interest. This brings about the possibility of 
calculating the interaction energy between vortices as a 
function of their distance. Then, by investigating confor¬ 
mations which minimize the corresponding energy, we are 
able to present a simple physical picture of the underly¬ 
ing phenomena which lead to elsewhere observed vortex 
states. 

This paper is organized as follows. In Sec. II, we briefly 
present the main aspects of the theoretical approach, fo¬ 
cusing on the Euler-Lagrange equations. Sec. Ill sum¬ 
mons our numerical results and displays a discussion con¬ 
cerning experimental realization in each considered case. 
While section III. A focuses on two-component BECs with 
contact interaction, section III.B concerns bound states 
of vortices in two- and three-component condensates fea¬ 
turing also a coherent coupling of the Rabi type. Our 
concluding remarks are presented in Sec. IV. 


II. THEORETICAL APPROACH 

We begin our consideration from the Gross-Pitaevskii 
energy functional for an Ac-component homogeneous 
BEG. When set into rotation above a given critical fre¬ 
quency, superfluids acquire angular momentum in the 
form of a vortex [301 El]. In turn, these can be repre¬ 
sented by nodes in the wave-function 'k Elj- Here we 
adopt a different approach and include vortices directly 
in the form of a node, disregarding rotation. There¬ 
fore, only contributions from the kinetic and interaction 
energies should enter in the energy functional. More¬ 
over, since vortices will be placed manually in fixed po¬ 
sitions, we also do not need to account for the energy 
provided by the angular momentum of the system, L^, 
which is only necessary to allow the vortex formation in 
the Gross-Pitaevskii theory for a rotating BEG. Then, in 
the present case, the energy density functional reads 

Nq f fc2 '1 

^ = + + (I) 

Q;=l ^ 

where the components are labeled by the index a. Here, 
the first and second terms, respectively, account for the 
usual kinetic and contact interaction energies within each 
component. In addition, the term Vij stands for the inter¬ 
component coupling energy density. 

In what follows, we will be interested in cases where 
Vij can be written as 

£5bI" - > (2) 

i j>i 


where the first term stands for the contact interaction 
and the last one represents the Rabi term. The lat¬ 
ter characterizes the internal coherent coupling between 
the components and has been already experimentally im¬ 
plemented by inducing an external driving field, which 
allows particles to move from one hyperfine spin state 
to another ESj- These interaction terms are controlled 
by the parameters ga, gij and Wij, respectively. How¬ 
ever, it is convenient to redefine them in order to have 
dimensionless units of energy and length. Therefore, 
for a two-dimensional system, we introduce the units 
Si = h?p\/2mi for energy density and ^ for 

distances. Here, pi is the average particle density of the 
first species. With E = SjSi and r = P'/^, the energy 
density can be written in its dimensionless form as 

+ - Wij (V'iV’i + V'iV'D] , 

i j>i 

where ipa. = is the dimensionless order parameter. 
This procedure leads to the definition of the mass ra¬ 
tio Mia = iTii/ma, the dimensionless contact interac¬ 
tion strengths "fa = da Pi/Si and jij — gijp\/Si, and 
the dimensionless Rabi frequencies LOij = w^pi/Si. The 
contact parameters can be well controlled in experiments 
by means of Feshbach resonances and appropriate values 
should enable the visualization of properties of interest. 

In order to measure the interaction potential between 
two vortices, we consider the total energy as a func¬ 
tion of the inter-vortex distance. This is justified be¬ 
cause the system is free from external contributions, so 
that the vortex position only affects the vortex-vortex 
interaction. Moreover, calculating the system energy 
as a function of the distance between vortex pairs re¬ 
quires the vortex fixing, which is achieved by fixing the 
27r whirl of the phase of the condensate wave function, 
thereby introducing a node in it. For one vortex in the 
a-component, for example, the wave function is given as 
'>Pa{x,y) = faix, y), with Na the number of 

particles in the a-component and the corresponding 
winding number, i.e., the quanta of circulation carried by 
the vortex. In addition, the angle 9k is defined around the 
locus of the fc-th vortex. Correspondingly, for two vor¬ 
tices, we have i/’a{x,y) = y). 

Since such vortex dimer structure is not circularly sym¬ 
metric, Cartesian coordinates are an appropriate choice, 
with is written as 

( I ■ \ nfc,c«/2 

Xk,a + iyk,a \ 

^k,Oi ^Uk^a J 

where rk,a = (a^fe.aj 2 //c.a,0) stands for the in-plane po¬ 
sition vector with origin at the center of the vortex k. 
As a consequence of the fixed circular phase around the 
vortex, singularities in the amplitudes appear naturally 
from the energy functional minimization. 
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We remark that the present ansatz is general enough 
to allow for considering components whose vortices might 
differ both in position and winding number. This is im¬ 
portant in order to obtain the total energy of the sys¬ 
tem in the presence of the inter-component coupling as a 
function of the relative distance of the vortices and then 
verify whether one has found the lowest energy configu¬ 
ration or not. Indeed, in the absence of inter-component 
coupling = ojij = 0, we have found that the single¬ 
vortex state is infinitely degenerate and the respective 
position of vortices in different components is irrelevant. 
This supports our statement above that the total energy 
depends only on the inter-vortex distance. In the pres¬ 
ence of an inter-component coupling, however, vortices 
in different components will rather stay on top of each 
other (separate away) for attractive (repulsive) effective 
vortex-vortex interaction potentials. 

With the adequate ansatz for the condensate wave 
function at hand, we minimize the total energy, which 
is obtained by integrating in space the quantity E — 
with respect to /„. Here, fia stands for the 
chemical potential of the the a-component and is intro¬ 
duced to keep the corresponding number of particles con¬ 
stant, which means that we are not taking population 
transfer between different hyperfine spin states into ac¬ 
count. The minimization process is implemented by nu¬ 
merically solving the set of W Euler-Lagrange equations 

^laEafa T ^afa fafct 

Afc Wc 

+EE (^fi^a.,j “t“ fj^a.,i') (5) 

i j>i 

-Nlij {fjfzSi^a + = 0 , 

where 


III. NUMERICAL RESULTS 

The numerical solution of the CGP equations is ob¬ 
tained considering a two-dimensional rectangular system 
divided in a uniform square grid 4000 x 2000 with to¬ 
tal dimension 1000^ x 500^, by using the finite-difference 
technique and a relaxation method suitable for non-linear 
differential equations. This leads us to the lowest-energy 
vortex structure that satisfies the constraint that vortices 
are placed in the fixed positions. The obtained conden¬ 
sate amplitudes are then substituted back in the energy 
density, which is numerically integrated, yielding the to¬ 
tal energy U = J E(^a)dxdy (dimensionless in units of 
El = fi^^) of the corresponding vortex configuration. 
With the total energy as a function of the vortex-vortex 
distance at hand, our theoretical approach will then be 
used to investigate vortex-state solutions of the Gross- 
Pitaevskii equation in the presence of rotation, reflecting 
the real experimental conditions under which vortices are 
generated. 

In what follows, we present the results first for BEGs 
with contact interaction only, where we focus on the two- 
component case, and then also for coherently coupled 
two- and three-component BEGs. For simplicity, we took 
the same number of particles per component for all situ¬ 
ations, N = 5 X 10"^, as well as the intra-component cou¬ 
plings (71 = 72 = 73 = !)■ Since the investigation of the 
vortex-vortex interaction is performed within the misci¬ 
bility condition: 7 i 72 ~ 7 i 2 > 0, the used inter-component 
coupling parameters are within the upper and lower lim¬ 
its I 712 I < 1. We denote by the total energy for the 
case of i vortices placed in the first component, j vortices 
placed in the second, and k vortices placed in the third 
(if considered). 


Qij = cos [inii9ii - nij 6 »ij) -h (n 2 i 6 » 2 i - n 2 j 02 j)]. ( 6 ) 


A. Two-component BEG with contact interaction 


At this point, we have chosen the same number of par¬ 
ticles for all components by setting Na = N. Notice 
that Eq. ([^ has the same form of the Gross-Pitaevskii 
equation and we shall call it constrained GP equation 
(GGP). One important difference, however, should be 
emphasized: the term 

(7) 


which arises from the kinetic contribution of the energy 
density functional, carries two extra terms and 
allowing us to manually set the coordinates of fixed vor¬ 
tices as well as their winding numbers. These terms are 
given by 


= 


^2,ct^2,a 
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Finally, we remark that a similar procedure has al¬ 
ready been sucessfully applied to superconductors in the 
framework of the Ginzburg-Landau theory [5^ . 


Let us first consider the case of a two-component con¬ 
densate, with two vortices placed at {x, y) = {±d/2, 0 ) so 
that the distance between them is d. We shall consider 
components with different masses, as this has been found 
to lead to unusual lattice conformations in Ref. m- By 
assuming a mass ratio of M 12 = 2.0 between the different 
components, we plot in Fig. [^the vortex-vortex interac¬ 
tion for different values of the inter-component coupling 
strength: 712 = —0.20 (circles); 712 = —0.45 (squares); 
712 = —0.75 ^iangles); and 712 = —0.98 (downward tri¬ 
angles). Fig. [^(a) and (b) represent the intra-component 
vortex-vortex interaction for the first and second compo¬ 
nent respectively, whereas Fig. [^c) depicts the inter¬ 
component vortex-vortex interaction. 

It turns out that for a two-vortex configuration, the 
interaction potential is always a monotonic function of 
the distance, where the attractive or the repulsive be¬ 
havior is determined by the properties of the particle- 
particle interaction. Nonetheless, one interesting fea¬ 
ture can already be identified at this level: for inter- 
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FIG. 1 . (color online) Vortex-vortex interaction potentials, 
considering M12 = 2.0, for (a-b) two vortices in the same com¬ 
ponent and (c) vortices in different components. Colours and 
symbols denote different inter-component coupling strength: 
712 = — 0.20 (blue circles); 712 = — 0.45 (yellow squares); 
712 = — 0.75 (red triangles); and 712 = — 0.98 (green down¬ 
ward triangles). 


component attraction, increasing the corresponding in¬ 
teraction strength I 712 I weakens the intra-component 
vortex-vortex repulsion. This can be understood in terms 
of the form of the amplitudes f{x,y). Indeed, as the at¬ 
tractive inter-component interaction becomes stronger, 
the energy is lowest for the largest possible superposition 
of the amplitudes, leading to a decrease in the density 
of one component in the spatial region where the other 
one has a vortex. Consequently, as I 712 I increases, these 
depletions should become more akin to vortices and the 
repulsive nature of the intra-component vortex-vortex in¬ 
teraction is softened by the depletion-vortex attractive 
interaction. 

Before we explore the consequences of the depletion 
induced in one component due to the presence of vortices 
in the other one, let us turn our attention to the mass 
ratio Mi 2 . According to the Feynman relation of the 
vortex density in a rotating superliuid [32], the vortex 
density is proportional to the particle mass, according to 


'^a^R 

^V.OL - Z 5 

TTfl 


( 8 ) 


where flu is the angular rotation frequency of the conden¬ 
sate. Thus, for a mass ratio of M 12 = 2.0, it is natural to 
investigate situations where two vortices (one vortex) are 
placed in the heavier (lighter) component. To this end, 
we have pinned a vortex in the second component in the 
center of the mesh and have placed two further vortices 
in the first component, separated by a distance d, as il¬ 
lustrated in Fig. [^by a contour plot of both components. 



FIG. 2 . (color online) The occupation number density contour 
plots for: (a) first component with two vortices positioned at 
(—d/2,0) and (d/2,0) and (b) second component with a sin¬ 
gle vortex in the center of the mesh. Notice the depleted 
density in one component in the positions where the other 
component features a vortex. The investigation of the sys¬ 
tem energy with respect to d for this conformation enables us 
to identify possible bound states in a two-component Bose- 
Einstein condensate with mass ratio M12 = 2 . 0 . 


In Fig. colors represent the density occupation number 
for both order parameters, and vortices coordinates are 
explicited as a function of the inter-vortex distance d. In 
this case, the system energy is no longer a vortex-vortex 
interaction potential. However, this conformation should 
enable the identification of possible bound states pro¬ 
vided by the competition between the inter-component 
and intra-component interactions, such as dimers and gi¬ 
ant vortices I32H1II in the same component. Notice the 
weaker density in one component in the positions where 
the other component features a vortex, characterizing the 
depletion effect discussed before. 

The total energy of this three-vortex configuration is 
plotted in Fig. |^(a) and (b), whereas Fig. [^c-f) demon¬ 
strates the previously discussed density depletion for sev¬ 
eral values of the relative interaction strength. 

From Fig. da), one sees that for low values of I 712 I 
the energy is a monotonically decreasing function with 
the distance d, indicating repulsion between the vortices 
in the same component. However, this behavior is com¬ 
pletely changed for 712 = —0.98, where the total energy 
turns into a monotonically increasing function, allowing 
vortices of the same component to occupy the same po¬ 
sition, thus, forming a giant vortex. Indeed, as one can 
see in the magnification of this result in Fig. [^b), the 
energy minimum corresponds to a vanishing distance be¬ 
tween the vortices in the same component. At this point, 
it is evident that the depletion effect is the main cause 
for the behavioral change of the effective interaction. In 
fact, the strong coupling between the components makes 
the intra-component repulsive contribution less relevant 
than the inter-component vortex-vortex interaction, and 
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FIG. 3 . (color online) (a) Total energy for three vortices con¬ 
figuration U21, characterizing a binary system with M12 = 2. 
The vortex in the second component is hxed in the center of 
the mesh, whereas vortices in the first component are placed 
at symmetric positions separated by distance d. (b) A mag¬ 
nification of the particular case 7 i 2 =- 0.98 shows that in the 
favoured conhguration, the two vortices in the first component 
are located on top of each other, and on top of the vortex in 
the second component, (c-f) Occupation number density of 
both components, considering vortices separated by 100 
for all cases considered in (a), where the solid (dashed) line 
stands for the first (second) component. 


the attractive effective interaction is achieved for every 
distance d. In Fig. [^c-f), the occupation number den¬ 
sity for vortices separated by 100 ^ shows the depletion 
effect, where depletions become more akin to vortices for 
higher values of inter-component particle coupling I712I. 
It is noticed that vortices become larger when 712 is in¬ 
creased, causing the vortex strong overlap illustrated in 
Fig if). This effect depends on the inter-component 
particle coupling and vortex density. 

Before going into further details concerning vortex 
core, let us first analyze how the effective interaction be¬ 
tween vortices changes from repulsive to attractive, by 
considering intermediate values of 712 between 712 = 
— 0.75 and 712 = — 0 . 98 . We have found that, for 
712 = — 0.90 and 712 = — 0 . 92 , as illustrated in Fig. ia), 
the energy curves still have repulsive characters, but they 
exhibit shoulders at finite distances, which could indicate 
the formation of non-triangular lattices in these minima 
of the interaction potential, despite the overall repulsive 
behavior. In Fig. ib), by setting 712 = — 0 . 95 , we have 
obtained a total energy which exhibits a non-monotonic 
behavior and acquires a clear minimum at finite inter¬ 
vortex distance. This strongly indicates that the presence 
of a vortex in one component can cause the formation of 
a bound state of two vortices (dimer) in the heavier com¬ 
ponent, illustrated in Fig. [^c). This results corroborates 
the states obtained in numerical experiments of Ref. m, 



FIG. 4 . (color online) (a) Total energy of a three vortex struc¬ 
ture in a two-component BEG, for 712 = — 0.90 (blue circles) 
and — 0.92 (red squares). The curve bumps may lead to non- 
triangular lattices, despite the overall repulsive behavior, (b) 
Particular case of a triple-vortex structure with 712 = — 0 . 95 . 
The minimum far from the origin allows for the formation of 
dimers for specific vortex densities, (c) Dimer configuration 
profile for 712 = — 0 . 95 . 


where the Abrikosov lattice of vortex dimers in one com¬ 
ponent, sitting on single vortices in the other, was found 
in a binary mixture. 

In addition to effects caused in the vortex-vortex in¬ 
teractions, leading to different lattice conformations, the 
depletion effect also leads to formation of vortices with 
different core sizes according to the choice of the 712 pa¬ 
rameter. Actually, this was pointed out before in Fig. 
I^c-f), where more negative values of 712 result in vor¬ 
tices with larger core sizes. To investigate this depen¬ 
dence, in Fig. by setting a vortex in only one of the 
two components, the vortex core radius was measured 
for different values of inter-component coupling, between 
712 = — 0.90 and 712 = 0 . 90 . This was done by measur¬ 
ing the distance between the center of the vortex and the 
point at which the order parameter decreases by half of 
its long-range convergence value. 

Since the components have different mass, in Fig. 
we investigated this effect for vortices in both compo¬ 
nents, where the core size for a first (second)-component 
vortex is presented by the solid (dashed) line. It is ob¬ 
served that the more the two components are coupled, 
the larger the vortex radii will be. This effect actually 
comes from the depletion originated by the vortex deple¬ 
tion which contributes to an increased vortex size. For 
both components, this behavior is qualitatively the same, 
suggesting that the different mass causes simply a shift 
between two curves. 

Let us at this point address the question of how re¬ 
alistic the present predictions are from an experimental 
point of view. For example, a two species BEG with 
a tunable interspecies interaction has already been pro¬ 
duced from the mixture of ®^Rb and for which one 
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FIG. 5. (color online) Vortex core size as a function of the 
inter-component coupling, for M 12 = 2.0. The blue solid line 
(red dashed line) represents the core radius of a vortex in the 
first (second) component. 

has Ml 2 « 2.1 [15]. In that particular study, a Feshbach 
resonance around a magnetic field of S « 79G allows for 
tuning 012 from positive to negative values. For example, 
for 80.7G 012 ~ — 185oo is reported, with oq the Bohr ra¬ 
dius. This clearly opens up the possibility of obtaining 
the appropriate negative values of 712 needed to exper¬ 
imentally observe the predicted vortex states. We also 
remark that, for a ®^Rb-®®Rb mixture, tunability of the 
interspecies interaction has already been used to probe 
various mean-field regimes such as spatial separation as 
well as the formation of long lived droplets (see Ref. [48] ). 

B. Vortex dimers and trimers in coherently 
coupled BECs 

In this subsection, we consider the presence of an inter¬ 
nal coherent coupling of the Rabi type between the differ¬ 
ent components. In order to do so, an extra term should 
be added to the energy density functional, as pointed out 
in Sec. |TT| For a system with Nc components, the Rabi 
contribution in Eq. ([^ can be rewritten in the form 

Na 

- 2^ ^ ujijl-ipiW-ipjl cos( 6 », -9j). 

i j>i 

If the signs of coefficients are all positive, the ground 
state energy is achieved when all the phases 0 i are the 
same. Recalling that the phases fix the position of the 
vortices, this means that vortices of different components 
would overlap, featuring an attractive potential. On the 
other hand, by setting positive values for the contact cou¬ 
pling 7 y , a repulsive inter-component vortex-vortex in¬ 
teraction emerges. The balance between these two types 
of interaction can lead to molecular conformations of vor¬ 
tices. Unlike the previous case, these bound states do 
not arise due to the mass difference between the compo¬ 
nents, but due to the competing interactions induced by 
the Rabi contribution. To illustrate the above argument. 



d/^ 


FIG. 6. (color online) Two-vortex inter-component interac¬ 
tion potential in the case of competing contact and Rabi cou¬ 
pling, cj = 2.1 X 10~® and M 12 = 1, for several values of 
712. 


consider, for instance, the case with two components with 
equal masses Mi^ = 1 , with one vortex in each compo¬ 
nent. 

As illustrated in Fig. a vortex dimer may be formed 
in the case of a short range repulsive and long range 
atractive inter-component potential, as observed for pos¬ 
itive values of "fij , where a minimum of interaction energy 
is found at finite inter-vortex distance due to a compe¬ 
tition between two types of inter-component interaction. 
On the other hand, negative values of 712 lead to an ex¬ 
tended attractive interaction in the long range. 

Let us now address the question of the formation of 
vortex trimers in a three-component condensate, i.e., 
bound states consisting of three vortices. In an at¬ 
tempt to reduce the excessive number of parameters, 
we shall consider the case with a;i 2 = W 13 = 0^23 = w 
and 7 i 2 = 713 = 723 = 7 . This choice also helps in 
finding the minimum energy configuration, since in this 
case, for symmetry reasons, it is expected that the three 
vortices are arranged equidistantly. Thus, we manually 
placed the vortices at the vertices of an equilateral tri¬ 
angle of side d, whose loci then read S*! = ( 0 ,-s/Sd/d), 
S '2 = (—d/ 2 , —i/Sd/d) and S 3 = (d/ 2 , —i/Sd/d). This is 
illustrated in Fig. where contour plots of the densities 
of the three components are shown. Notice the presence 
of the depletion effect in this case as well. Subsequently, 
we investigated the minimum energy conformation as a 
function of d for several values of w. 

In Fig. I^a), where we considered 7 = 0.45, there is 
indeed a finite value dmin of the triangle side that mini¬ 
mizes the interaction energy Um. Moreover, Fig. §(b) 
and (c) display power law dependencies of dmin in w, for 
7 = 0.45 and 7 = 0.60, respectively, which explains the 
vortex trimer configurations observed e.g. in Fig. 3 of 
Ref. |20|. Despite the fact that the exponents are quite 
similar in the two considered cases, they should depend 
on both 7 and number of particles N. Indeed, consid¬ 
ering the extreme case where 7 vanishes, the inflection 
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FIG. 7. (color online) The occupation number density con¬ 
tour plots of: (a) first component with a vortex positioned at 
(0, \/3d/4); (b) second component with a vortex positioned 
at (—d/2, —y/id/A) and ; (c) third component with a vortex 
positioned at (d/2, —-\/3d/4). 



FIG. 8 . (color online) (a) Total energy [/m, for the equilat¬ 
eral triangle vortex configuration as a function of the side d for 
7 = 0.45 and for several values of uj. (b) Plot of the optimal 
side of the equilateral vortex triangle against uo for 7 = 0.45. 
The red line represents a power law with exponent —0.3214 
and coefficient 0.9265. (c) Same as (b) but for 7 = 0.60. The 
red line represents a power law with exponent —0.3186 and 
coefficient 1.0746. 


point of the vortex interaction should always be at d = 0 
for any value of w, leading to a vanishing exponent in the 
power law dependence dmin (7)- 

Since the geometry was imposed from the beginning 
by setting the vortices at the vertices of the equilateral 
triangle, the existence of a value dmin which minimizes 
the total energy does not guarantee that this configu¬ 
ration corresponds to the ground state. In order to al¬ 
low for different configurations and thereby check if the 
equilateral triangle is really a minimum of the energy, 
we have allowed for different positions of the vertex ^i, 
while keeping S2 and S3 fixed for w = 1.7 x 10 “"^ and 
7 = 0 . 45 . As can be seen from Fig. the equilateral tri¬ 
angle configuration turns out to be stable, corresponding 
to at least a local minimum of the interaction energy. 

Although not very common, there are other ways to 
obtain bound vortex states which do not arise from a 
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FIG. 9. (color online) (a) Total energy Fm as a function 
of the displacement of the first component vortex Xvi from 
5i position in the x direction, (b) Total energy Fin as a 
function of the displacement of the first component vortex 
from the origin in the y direction. For both cases, we have 
considered 7 = 0.45 and u — 1.7 x lO”'^. The vertical dashed 
line represents the coordinates of first component vortex for 
the minimum energy equilateral triangle configuration. 


competition between different types of interaction. In 
fact, it is also possible to obtain dimers and trimers 
exclusively from Rabi coupling in condensates with at 
least three components. These vortex states arise from 
frustration between the phase locking tendencies. As 
it was mentioned previously, a set of positive Rabi fre¬ 
quencies W12 = Wi3 = W23 = w > 0 would lead to 
di = O2 = ^3 as a minimum of the energy, causing the 
same spatial occupation for all three vortices. However, 
by assuming, for example, a;i2 = 0J13 = —1023 = w, it is 
not possible to satisfy the minimization energy condition 
cos((()i - cj)2) = cos((()i - (/13) = I and cos((/)2 - 03 ) = -1 
simultaneously, which characterizes a frustrated system 

mm- 

Since the frustration phenomenon arises basically from 
the problem of minimizing the total energy with respect 
to the order parameter phases, we must redefine our 
ansatz, adding extra phases to our order parameters in 
order to account for the phase tendencies which mini¬ 
mize the energy. Thus, by assuming an ansatz of the 
form ipnix.y) = •\/iVe*®"e“*'^"/„(a;,y), the consequences 
of choosing Rabi frequencies which lead to a frustrated 
system were investigated for a;i2 = W13 = —0:23 = w 
and 712 = 713 = 723 = 0, where we have considered 
01 = 02 = 0 and several values for 03. 

In Fig. [Iota), we pinned vortices of the first and sec¬ 
ond components on top of each other, in the center of the 
mesh and displaced the vortex of the third component at 
distance d from the origin. This procedure was made for 
different values of phase 03 . The obtained results show 
a potential that suggests the formation of a bound vor¬ 
tex state, with the third component vortex outside the 
origin, for 03 = 0, where the potential minimum is at 
d « 4 ^. This conformation arises from the competition 
between the attraction with the vortex of the first com¬ 
ponent and the repulsion with the vortex of the second 



















FIG. 10 . (color online) (a) Total energy Uwi, with vortices of 
the first and second components pinned in the center of the 
mesh, as a function of the third component vortex distance 
d from the origin, for different values of the third component 
extra phase 0 (solid blue line), vr/O (green dashed dotted 
line), 7r/3 (yellow dashed line) and 7r/2 (red dotted line), with 
a>i2 = Wi3 = — aj23 = 1 X 10 “^. (b) Total energy ? 7 iii for 
the same conformation of (a), keeping <()3 = 0 and assuming 
different values of ta: cu = 1 x 10“^ (blue circles), uj = 2x 10“^ 
(green squares), cu = 3 x 10 “^ (yellow upward triangles) and 
ui = 4 X 10 ~^ (red downward triangles). 


component. By symmetry, the same result should also 
appear for (j) = tt. This configuration has the lowest 
energy among the other investigated cases, however, we 
can not guarantee that this conformation is the ground 
state. Actually, finding the ground state conformation 
would require the variation of all three extra phases and 
also other vortices positions in the grid, which is very 
expensive computationally and is left for future investi¬ 
gation. Nevertheless, these results are sufficient to ensure 
that vortex molecules are stable states. In Fig. [Mb) , we 
show that increasing the Rabi coupling makes the frus¬ 
tration effect even more robust and leads to a deeper 
minimum in the interaction potential. 

In view of the present predictions concerning Rabi cou¬ 
pled BECs, it becomes important to consider the pos¬ 
sibilities of experimental realization of adequate sam¬ 
ples. In that respect, we focus on the phenomena of 
interest, namely Rabi coupling and homogeneous con¬ 
finement. They have both already been experimentally 
achieved with ®^Rb. Before we turn to the particular ex¬ 
perimental setups, we recall that in our calculations the 
parameter governing the Rabi coupling is w = wpijSi, 
where w = hfl relates the Rabi energy w to the Rabi 
oscillation frequency fl. 


We first consider a recent study in Ref. [H], demon¬ 
strating long-range coherence in quasi-two-dimensional 
Bose gases. There, two-dimensional densities of about 
a few hundred ®^Rb-atoms per square micrometer are 
achieved and combined with harmonic trapping in the 
third direction whose frequency Vz ranges from 350 up 
to 1500 Hz. Moreover, for ®^Rb, the (three-dimensional) 
s-wave scattering length is Ug « IOOuq, with oq the Bohr 
radius. For the values of the Rabi couplings we turn to 
Ref. [ 33 ], where Rabi oscillation frequencies of the order 
27 r X 10 ^5“^ have been realized. We estimate the Rabi 
coupling parameter w by calculating Ei = for these 
values of the particle density and Rabi frequencies, to ar¬ 
rive at UJ = = 2 ^ ~ 1.7 X 10 “^, where we have 

associated a length a^iabi = with the Rabi os¬ 

cillation. This value is already higher than ones needed 
for the effects discussed in this manuscript. 

It is also possible to improve the estimates of the actual 
experimental values of the coherent coupling by including 
the effect of the transversal harmonic trapping, render¬ 
ing the calculation more realistic. Indeed, in order to 
obtain typical values of w in quasi-two-dimensional sys¬ 
tems, one should correct the s-wave scattering parameter 
g for the freezing of the third dimension. This is done 
by dividing the three-dimensional gsD = 47rh^as/mi by 
\/^az, with Oz being the oscillator length in the third 
dimension [JS]. If one then associates the energy Ei 
with the contact interaction energy gpi, a slightly dif¬ 
ferent expression for the coherence parameter emerges 
UJ = 2 ^ — « 1 X 10 “^. Notice that, besides 

“RabiPl VStTOs ’ 

the fact that this value of uj is even higher in quasi-two- 
dimensional traps, this expression allows to identify the 
trap frequency in the third dimension, in addition to the 
Rabi oscillation frequency and the particle density, as 
tuning knobs for studying the effects of coherent coupling 
in multi-component BECs. 


IV. CONCLUSION 


In summary, we have presented a method for calculat¬ 
ing pairwise vortex-vortex interaction potentials in multi- 
component Bose-Einstein condensates, capable of reveal¬ 
ing the underlying reasons for the unusual vortex con¬ 
figurations, such as different lattice geometries and few- 
vortex clusters. We have applied our theory in specific 
examples and could thereby clarify the formation of a 
lattice of dimers in rotating two-component BECs as 
well as the dimer/trimer state in a coherently coupled 
two/three-component BEC by pointing out the role of 
the non-monotonic interaction potential with respect to 
the inter-vortex distance in both cases. Our analysis of 
the present day experimental capabilities, combined with 
simple estimates of relevant parameters, indicates that 
considered effects are on the verge of being directly ob¬ 
served in the lab. We remark that the present theory can 
be straightforwardly adapted to any number of conden- 
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sate components, providing a tool to anticipate different 
aspects of vortex physics in Bose-Einstein condensates. 
As a matter of fact, even for one-component BECs fea¬ 
turing the anisotropic and long-range dipolar interaction, 
vortex phases such as bubbles and stripes have been fore¬ 
seen [ini HH |49] , for which a thorough understanding of 
the role of vortex-vortex interaction is still lacking. 
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